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ABSTRACT 

A  m«  aoUtioa  is  pwnlii  far  tko  maaUmlloa  of  proilaa  o f  aareaol  vokm  oliittiaa  coaffi- 
ckau  frawi  tka  aofcy  fcathtiKiwf  rotanaa  af  a  aaoaaatatk  Ma«W-war%laa|«k  lidar  tyaUm.  TUa  iavaraa 
pwblaaa  k  aohrad  by  wtiHsiaf  a*  lafonaatka  tkaeaalk  aaaUiad  baaad  am  Um  priacipia  af  amiwaaa  ctom- 
aatropy  (MCE),  wkkk  npcaaaala  aa  abjartiva  aad  ntiaaai  approach  for  Um  affactiva  iacorporatioa,  into 
Um  bvanioa  procadaio,  of  both  prior  iafaraaatiaa  ia  tka  form  of  aa  iaitial  aatiaaata  of  tka  axtiactioa 
coaflkiaat  ud  additional  bformatfea  M  Ika  lira  af  Ika  okatd  lidar  data.  A  akapia  aad  rokaal  a*, 
markal  ptwcadara,  kaaad  oa  Um  afltpaoad  algorithat,  la  davaiopad  to  compote  Um  MCE  raceaalraciioa 
af  Um  antiactioa  faactioa.  A  itak.'  af  aaaarkal  aaaiplaa,  baaad  aa  aoiay  rymtkatk  Bdar  data,  ara 
amployad  to  daaioaatrata  aad  avalaate  tka  atilky  aad  aftcary  of  tka  invaraioa  nMtkod. 
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INTRODUCTION 

i.  Aa  optk  I  nmota  i— h|  lyKia  bawd  oa  tka  light  dstertiag  tad  ranging  (lidar)  principle  can 
ba  atiliaed  for  tka  tat  car  ary  of  mini  properties  (a.g.,  mam  sad  ■■■bar  rnariatrttina.  partkle  tiac 
distribntioa,  aad  optical  coaataata)  of  aaroaol  doada  from  tka  backacattand  rat  arm  of  a  coatrolled  laser 
palaa.  Siace  praaaat  Bdar  ayatawa,  rack  aa  tka  Laaar  Cload  Mappar  (LCM)  (1),  an  cap aMa  of  |aaarating 
vary  akort  laaar  pa  lata  at  rapid  rapaUtioa  rates,  tkaaa  ayataaaa  effectively  poaaaaa  tka  capability  for  rapid 
aad  datailad  probiag  aad  r  weillaacs  of  aaroaol  cloada  wRk  aakaacad  taatithrity  aad  raaohatioo.  I  ad  aad, 
tka  com  pact  aaaa,  mobilty  aad  apaad  of  oparatlaa  aflordad  by  tka  LCM  araksa  it  a  nhubk  tool  for  tka 
rapid  monitoring,  ckaractarisatlea  aad  aaaaaaaaaat  of  tkraat  aaroaola  wkfck,  aa  a  coaaeqaaaee,  caa  provide 
tka  aaar  witk  a  mack  eartisr  warning  of  potaatiaBy  kaaardoaa  aitaatioaa  tkaa  caa  ba  provided  by  local 
moaitoriaf  procadaraa  baaad  oa  varioaa  diract  ia  aita  aanpUaf  matkoda.  Bowaver,  it  Ja  -aa  port  amt  to  not* 
that  a  diract  (local)  aaaaiag  tacaaiqaa  provides  a  aon  comprakaaarw  characterisation  of  tka  phyakal 
propart  iaa  of  tka  aaroaola  since  it  ia  baaad  oa  tka  actaal  c  apt  art  aad  aaamlaatioa  of  tka  partkalataa, 
wkereas  aa  optkal  ramota  sensing  tackaiqaa  maat  iaftr  tkaaa  propartka  from  tka  backacattand  pattanu 
of  tka  eiactromagaetk  waves. 

2.  Iadaad,  tract  atawapkaric  partkalataa  iataraet  witk  optical  radiatioa  ia  a  compikatad  maantr, 
tka  problam  of  inferring  phyakal  prcyerties  of  atmoapkark  aamaok  from  backacattand  radiatioa  ia  not 
straightforward  aad  ia  fraagkt  witk  aamaroaa  diflcsRiaa.  Coaaaqaaatly,  than  still  doss  not  exist  a 
totally  satisfactory  solatia*  to  this  problem,  despite  tka  eaaaidarehlt  rasaarck  effort  expended  on  it 
otar  tka  past  thirty  years,  Tka  di Realty  ia  formalatiag  aa  appropriate  sohitioa  to  tka  lidar  inverse 
problem  items  from  tka  fact  that  tka  rseoaatractioa  of  tka  optkal  proptrtiaa  of  aa  aaroaol  cload  from  a 
limited  sat  of  noisy  Udar  ratara  data  k  aa  i abate* Uy  iO-poaad  problem  ta  tka  tease  tkat  tka  aolation  is 
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Mk-nufM  (vis.,  tWt  ars  infinitely  many  wtkrtwi  fnetioM  tkat  an  cons tot— >  with  Um  |ma  data) 
aad  unstable  (rk.,  th*  aolatioa  do*  not  depend  tanrtananaly  on  the  data).  Consequently,  tka  formal 
malhaanatfcal  solution  for  tka  lidar  iiana  problem  [2}  vaa  fonnd  ta  m  plafond  by  bwtabilitfoa  aad 
ia  accuracies  wkkk  affectively  rendered  sack  a  eolation  aaalan  for  practical  applications  la  fact,  tka 
failure  of  tkia  exact  (La.,  foraul)  aolatioa  waa  aot  well  aadaratood  aad  wan  originally  attribatad  to  tka 
physical  affacta  arising  from  tka  a »flact  of  tka  makipfo  acattariaf  caatribotioae  ia  tka  treatment  ratkar 
tkaa  from  tka  inherently  ill- posed  aatare  at  tka  prokhm. 

L  Tka  bdar  invars*  problem  caa  be  uniquely  cobrad  provided  oaa  fo  fh*m  aa  ideal  aad  complete 
■at  at  lidar  bocfcacatteriag  returns,  bat  practical  Haatatioao  impoaad  bp  tka  mamoramant  epetam  a* 
trail  m  aooaa  rwaapiaf  simplifications  ariaa|  from  tka  aoffoct  at  makipfo  acattariaf  coatribatioas,  tarra 
to  tkwart  tka  attainment  at  tka  exact  mathematical  aolatioa.  It  ekoald  be  empkaafoad  tkat  ander 
realistic  drcomstancas,  tka  exact  aolatioa  caaaot  be  mad  to  completely  ckaractariae  tka  trae  propartiaa 
of  tka  aaroaol  clouds  from  tka  obearrad  noisy  aad  dfocrata  bdar  data,  Indeed,  tka  exact  aolatioa  doea 
aot  area  make  aay  aUoaraacm  for  tka  aoim  ia  tka  data  ia  aap  optimal  manner.  Nevertheless,  all  tka 
aolatioa*  propomd  for  tka  iararaioa  of  tka  Bdar  aqaatioa  to  date  kora  baaa  baaad  apoa  aom*  appropriate 
modification  of  tka  exact  mathematical  aolatioa.  Kfott  [S]  akomad  tkat  tka  exact  bdar  ia  vena  aolntion 
caa  be  atabiliaad  if  oa*  ckooem  tka  boaadary  vahm  (La.,  tka  value  of  tka  axtiactioa  fa  action  at  aom* 
specified  raaga)  at  tka  far  aad  of  tka  aaroaol  dead  ratkar  tkaa  at  tka  aaar  aad.  However,  tkia  required 
tke  determiaatioa  of  tke  boaadary  value  wkkk  caaaot  be  obtained  from  the  data  directly.  Even*  [4] 
demonstrated  bow  tke  Utter  problem  caa  ba  owiioma  by  utibsng  aa  appropriate  calibration  of  tka 
bdar  system  fot  coajeactioa  with  a  aim  pi*  modification  of  Kfott’*  iavaraioa  procadora.  Finally,  it  ekoald 
bo  meatfaaad  tkat  two  popolar  bdar  iamtioa  procedure*,  aamaly  tka  slop*  aad  ratio  matkoda  (S),  caa 
ba  coaaidarad  to  ba  simply  first  order  apprexfoeaiioas  at  tka  exact  solutions. 

4.  Ia  lifkt  of  tka  nonaarons  difficulties  mcoanterad  ia  tka  formahiioa  of  aa  appropriate  aolatioa  for 
tka  lidar  inverse  problem,  it  fo  important  to  raaaamina  tka  problem  from  tka  point  at  view  at  evaluating 
tke  physical  validity  aad  meaaingfmfoeae  at  tka  recovered  aolatioa  ia  relation  to  tka  not  or*  aad  quality 
of  tka  data.  Aa  stated  previously,  although  formal  analytical  solutions  exist  for  tke  inversion  of  tke  lidar 
aquation,  tke  difficulties  that  arias  ia  tka  baplamaatatioa  at  them  eolations  reside  ia  tke  fact  that  the 
bnckscattered  radiation  is  aot  known  precisely  ewer  the  entire  rang*  interval.  ladaad,  tka  finite  and  noisy 
character  inherent  in  nil  realistic  lidar  data  ultimately  contributes  to  tka  ifl-po**d  nature  of  tka  lidar 
invars*  problem,  la  this  paper,  aa  information-theoretic  method  fo  developed  for  tke  inversion  of  the 
lidar  equation.  Tkia  method  pern  its  tke  objective  incorporation  of  prior  information  obtained  from  the 
physical  character  of  the  extinction  function  aad  of  additional  information  encoded  in  tka  lidar  data 
Tka  interaction  between  the**  two  component*  of  information  ia  eimpfo  :  prior  kaow ledge  relating  to 
varions  physical  characteristics  tkat  tke  aolntion  meet  pneama  fo  utilised  to  eliminate  a  vast  majority 
of  member*  from  tke  infinite  set  of  eolation#  consistent  with  the  finite  uncertain  data  and  from  tka 
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rein  awing  cwJiiUtM  in  thk  eolation  Mi,  ik«  principle  of  parsimony  k  etilissd  to  select  tka  aimpleat 
eolation,  Lt.  tka  solation  wkidk  posse  mm  tka  laaat  amount  of  detail  or  information  til  at  waa  not  already 
known  or  expected.  Tka  bane  philosophy  bakind  Making  tka  aimplaat  modal  (aohition)  raaidaa  in  the 
fact  tkat  one  doaa  not  want  to  ba  mialad  by  any  axtraaeoas  faatvraa  in  tka  aohition  that  are  not  eaaantial 
in  matching  tka  data.  Alternatively,  it  ia  axpactad  tkat  tka  feat  area  which  appear  in  tka  panunonoua 
eolation  are  significant  ia  tkat  they  an  (pacifically  required  by  tka  data.  In  tkie  paper,  tka  application 
of  tkma  idaae  to  tka  aointioa  of  tka  lidar  invane  problem  ia  quantified  in  tka  fern*  at  tka  principle  of 
mininmm  otom  entropy. 

THEORY 

S.  Tka  eingW-acat taring  lidar  eqaatioa  wkick  relatea  tka  voiame  extinction  coefficient  to  the 
backacattered  radiation  for  a  monoatatk,  nngla- wavelength  lidar  ayatam  ia  given  by 

where  r  ia  tka  range,  P(r)  ia  tka  backacattered  power,  Pb  ia  the  initial  laser  poke  power,  A  ia  tka  receiver 
area,  e  ia  tka  velocity  of  light,  r  ia  tka  la ear  poke  time  deration,  fi[r)  k  tka  bachecattering  coefficient 
function,  and  a,(r)  ia  the  volume  extinction  coefficient  function.  By  normaliaing  tka  lidar  ret  urn  power 
by  the  range  and  tka  ayatam  const  acta  to  form  tka  derived  aignal  S(r)  a  Ia[2raP(r)/(cr.PoJ4)],  it  ia 
poMible  to  rewrite  tka  lidar  equation  ia  a  form  which  removee  tka  effects  of  geometric  attenuation  and 
pure  ayatam  factors: 

S(r}~^(r)]-2jTa.(rVr'.  (la) 

A  La  actual  maaauremanta,  the  observed  backacattered  intensity  ia  corrupted  by  noke  and,  conse¬ 
quently,  the  measurement  model  is  taken  as 

IK  *  F(ri)  +  U,  »  »  1,2, (16) 

Hera  ja  ia  the  measured  logarithmic  ranga-aormalisad  backacattered  power  corresponding  to  the  range 
r,.  The  «,’•  are  measurement  and  process  error*  associated  with  the  jg’s  whose  standard  deviations 
are  denoted  as  <r<.  Tka  a,  which  characterise  the  uncertainty  in  tka  experimental  measurement  of 
the  jg’e,  are  assumed  to  be  known.  It  k  important  to  remark  that  these  errors  may  include  process 
contributions  arising  from  multi-scattering  effects  as  well  as  measurement  contributions  arising  from 
system  Mnsitrvities,  non-ideal  detector  respoaM,  and  numerical  and  digitisation  errors.  Note  that  eq. 
(la)  contains  two  unknown  functions,  namely  fi{r)  and  a«(r).  In  order  to  make  further  progress  on  the 
recovery  of  these  functions  from  the  observed  data,  it  k  necessary  to  make  an  assumption  concerning 
the  functional  relationship  between  the  bachacattering  aad  extinction  coefficients.  There  k  experimental 
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•ad  theoretical  evidence  [3,9]  that  under  a  widt  rmage  of  atmospheric  conditions,  these  two  coefficient 
functions  an  approximately  related  aa 


fit')  -  *(r)°»(r)>  (le) 

wksre  x(r)  ia  a  specified  fraction  o f  raaga  r  aad  the  power  k  far  most  aerosol  clooda  of  iataraat  ia 
coastrained  ia  the  interval  0.97  <  k  <  1.0.  Witk  tkia  assumed  faactioaal  dapaadaaca  between  fl(r)  and 
a.(r),  tka  lidar  ia pent  problem  involves  tka  recovery  of  tka  extinction  coefficient  from  a  given  aat  of 
noisy,  discrete  lidar  returns. 

7.  Now,  rappoaa  oaa  ia  givaa  or  kaa  available  by  aa  aharaativa  means,  soma  initial  (prior)  estimate 
p(r)  of  tka  actual,  but  ankaown  axtiaction  fn action  ct,(r)  (r  €  D  when  D  is  tka  support  of  the  extinction 
function,  vis.  the  set  of  all  poiats  of  range  r  wkara  a,(r)  is  non- varnishing).  Naturally,  this  initial  estimate 
is  selected  to  reflect  tka  user's  prior  information  on  tka  axtiaction  function.  This  prior  information  may 
coma  from  initial  modelling  stadiaa  of  generic  regional  aerosols  for  various  atmospheric  conditions  or 
from  tka  result  of  a  previous  inversion  obtained  using  a  different  data  set.  Now,  if  new  information  is 
subsequently  obtained  about  <a,(r)  in  the  form  of  backseat tered  return  data  from  some  lidar  experiment, 
then,  ia  light  of  tka  lidar  and  measurement  models  (et.  aqs.  (la),  (lb)  and  (lc))  which  relate  the 
extinction  coefficient  to  the  measured  data,  it  is  possible  to  specify  that  the  unknown  extinction  function 
must  ba  soma  element  of  aa  admissible  set  II  determined  aa  follows: 

n*  |g(r):5(r)=»ln[#e(r)]+*ln[g(r)]-2jr  g(r')A',  £(*  ~  S  X«}.  (2) 

where  x*/  is  some  prescribed  constant  that  reflects  the  level  of  noise  corrupting  the  data.  Recall  that  er, 
is  the  error  in  the  measurement  of  the  i-th  datum.  Observe  that  in  the  specification  of  the  admissible  set, 
the  effect  of  the  aoisv  observations  has  been  explicitly  accounted  for  by  the  usual  weighted  least-squares 
misfit  criterion 

imi 

which  ia  utilised  to  assets  the  goodaeaa-of-fit  of  the  theoretical  to  measured  backseat  tered  signals.  In  view 
of  uncertainties  in  the  measurements,  the  misfit  is  specified  to  be  leas  than  or  equal  to  some  predetermined 
constant  Xm-  Naturally,  if  the  observation  errors  constitute  an  uncorrelated,  sero-mean,  Gaussian 
process,  then  x3  poem  sees  a  clu-equare  distribution  aad  this  constant  can  be  chosen  a s  x«  *  Af,  the 
expected  value  of  x3.  Actually,  this  constraint  can  be  relaxed  somewhat  if  one  elects  to  choose  instead 
x\t  <  Af  +  l-JliA  (i.e.,  the  98  percent  confidence  limit  for  x3).  However,  it  is  important  to  emphasise 
that  a  statistical  model  for  the  observation  errors  is  not  required  for  the  ensuing  analysis.  Indeed, 
in  the  absence  of  aa,  a  priori  information  on  the  probabilistic  structure  of  the  measurement  errors, 
it  is  convenient  to  interpret  x3  simply  as  a  goodneaa-of-fit  criterion  that  measures  the  fidelity  of  the 
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model  date  to  the  observed  data.  Then,  the  selection  of  Xu  m  M  wonld  ba  interpreted  simply  aa  a 
root-meaa-sqaare  (SMS)  error  ol  1  ia  tba  fit  of  tba  modal  to  tba  obaanrad  data. 

8.  Having  specified  tba  praciaa  farm  of  tba  admbsibb  sat  II,  sappoat  mow  that  tba  iaitial  (prior) 
aotimata  p(r)  for  tba  extinction  function  ia  no  longer  consistent  with  tba  information  provided  by  the 
observed  data,  via.  p{r)  £  II.  Since  II  still  contains  infinitely  many  extinction  functions  that  are 
compatible  with  tba  obaanrad  lidar  data,  tba  problem  now  canters  on  tba  choice  of  tba  'beat'  estimate 
q{r)  of  ou(r)  from  II.  This  choice  should  not  only  ba  baaed  on  tba  given  data  as  embodied  ia  tba  form  of 
tba  admieeibla  sat  II,  but  should  also  incorporate  tba  properties  and  characteristics  of  tba  prior  estimate 
p(r)  in  soma  optimal  manner.  It  is  proposed  that  this  choke  should  ba  baaed  on  an  information-theoretic 
concept  usually  referred  to  as  the  principle  of  minimum  crose-entropy.  Whan  applied  to  the  present 
problem,  this  principle  dktataa  that  the  final  (posterior)  estimate  q(r)  for  the  extinction  coefficient 
should  be  cboaen  aa  that  member  of  the  admissible  set  II  that  poaseeees  the  least  arose -entropy  with 
respect  to  tbs  initial  (prior)  estimate  p(r).  In  otLer  weeds,  the  final  estimate  <?(r)  is  chosen  to  miniin  be 
the  cross-entropy  between  <?(r)  and  p(r)  defined  aa 

ff(«.P)“/  ?(r) la[9(r)/p(r)]dr, 

Jo 

subject  to  the  constraint  that  q  6  II  (vis.,  the  final  estimate  should  be  consistent  with  all  the  available 
information  contained  in  the  observed  data).  It  is  interesting  to  note  that  the  principle  of  minimum 
croee-entropy  b  referred  to  in  the  literature  by  a  plethora  of  names  depending  on  the  particular  field 
of  study  and  on  the  context  in  whkh  it  b  used,  including  such  terms  as  discrimination  information, 
1-divergence,  directed  divergence,  relative  entropy,  KeUback-Leibler  information  measure  and  expected 
weight  of  evidence  [7-12].  Furthermore,  it  should  be  noted  that  thb  principle  can  be  considered  to  be  the 
natural  generalisation  of  the  maximum  entropy  (MAXENT)  principle  which  was  originally  developed  by 
Maxwell,  Boltsmann  and  Gibbs  in  tbs  context  of  statbtkal  mechanics  and  later  extended  by  Jaynes  to 
the  status  of  a  general  inference  procedure  of  whkh  the  statistical  mechanics  application  was  merely  a 
special  case  [13-15).  In  particular,  the  principle  of  minimum  cross-entropy  reduces  to  Jaynes’  MAXENT 
procedure  if  the  prior  estimate  p(r)  b  selected  to  be  a  constant  or  uniform  function  over  the  domain  D. 

9.  The  rationale  for  the  minunnm  cross-entropy  principle  as  an  appropriate  criterion  of  choke  b 
elegantly  embodied  in  the  axiomatk  formulation  of  Shor*  ■-,d  ’  n  (11),  who  demonstrated  that  the 
principle  of  minimum  cress  entropy  b  ths  only  inference  procedure  thni  b  consistent  with  four  basic 
criteria  (axioms)  that  any  rational  method  of  inference  must  possess.  These  conabtency  conditions 
spscify,  among  othsr  things,  that  ths  solution  provided  by  any  rational  inference  procedure  must  be 
unique  and  invariant  under  coordinate  transformations.  From  another  viewpoint,  the  croee-entropy  can 
be  interpreted  as  a  measure  of  the  information  ‘distance*  between  the  initial  (prior)  and  final  (posterior) 
estimates  p(r)  and  q(r),  respectively.  Along  thb  vein,  a  geometrical  interpretation  of  the  minimum  cross- 
entropy  principk  can  be  depkted  as  ia  Fig.  1.  Thb  figure  shows  ths  final  estimate  q  for  the  extinction 
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function  aa  tl m  iifora*tk»-  or  I-prqjoction  o f  the  initial  mtimat*  f  onto  th*  admhribla  sat  H.  Rom  this 
ptnptctm,  Um  erwaantropy  B(q,p)  ua  ba  iitvpntU  ae  Um  ‘distant**  htmn  Um  U  >  estimate* 
p  ud  f.  Sine*  f  k  n  Um  ^projection  of  p  onto  IX,  this  distance  corresponds  to  th*  shortest 

distant*  between  p  aad  may  poiat  of  II,  vis. 

10.  Haaca,  B{q,p)  caa  ba  considered  to  ba  a  maaanra  of  th*  amount  of  information  encoded  in  th# 
data  tkat  k  not  already  ambaddad  ia  tba  initial  aatimata  p.  Indaad,  tua  agar*  indicate*  that  the  croaa 
entropy,  whan  viawad  aa  aa  information  distance  batwaaa  two  extinction  function*,  satkfie*  a  triangle 
equality  of  the  form 

B(a*,p)  *  B(a„q)  +  B{q,p), 

a  fact  that  caa  ba  raadily  verifed  by  aoma  limpla  manipulation*.  Whan  looked  at  from  thia  perspective, 
it  ia  readily  apparent  that  the  hal  aatimata  q  ia  doaar  to  the  tree  but  an  known  function  a«  than  ia  the 
initial  aatimata  p  since  H(a.,q)  <  H(a,,p).  Indeed,  the  final  estimate  q  is  the  aatimata  doeeat  to  the 
initial  aatimata  p  that  at  th*  same  time  verifies  the  data  constraint.  Finally,  it  should  ba  noted  that 
if  the  data  constraints  contain  no  new  information  on  the  extinction  function  other  than  that  already 
embedded  ia  th*  initial  estimate  p,  th*  final  estimate  q  ia  this  case  (aa  selected  by  th*  MCE  principle) 
ia  simply  th*  initial  eatimat*  p.  This  would  correspond  to  the  situation  where  p  is  a  member  of  the 
admiasibl*  set  II  aad  th*  I- projection  of  p  on  II  would  amply  b*  p. 

11.  Ia  summary  then,  whan  provided  with  aa  initial  (prior)  eatimat*  p(r)  of  th*  true  extinction 
coefficient  a,(r),  the  final  (posterior)  eatimat*  q(r)  is  chosen,  ia  adherence  to  th*  principle  of  minimum 
cross  entropy,  aa  that  eatimat*  which  solves  th*  following  functional  minimisation  problem  : 

*}  (3a) 

subject  to  th*  constraint 

?(*•)€  n,  (3») 

where  II  is  th*  admiasibl*  set  defined  in  *q.  (2).  It  is  important  to  note  that  since  II  possesses  the 

geometrical  structure  of  a  closed  aad  convex  set  ia  th*  muhi-dimensional  model  space,  there  exists  a 

unique  q{r)  €  II  (i.*.,  unique  up  to  a  set  of  measure  tero)  which  solves  the  minimisation  problem  posed  in 
eqs.  (3a)  and  (3b)  [10).  Additionally,  it  should  also  be  remarked  that  the  croas-eutropy  functional  H(q,p) 
is  positive  aad  convex  in  both  q  aad  e.  Observe  at  ,hk  point  that  the  reconstruction  of  th*  extinction 
coefficient  based  on  the  MCE  principle  requires  th*  solution  of  a  nonlinear  functional  (i.e.,  infinite* 
dimensional)  minimisation  problem.  Unfortunately,  there  does  not  exist  a  closed-form  solution  for  this 
problem  ud,  consequently,  some  numerical  method  must  b«  i,  cified  for  computing  u  approximation 
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to  tkk  solution.  The  development  of  a  numerical  algorithm  to  solve  aa  appropriate  iaite-dimanrional 
approximation  to  tkk  problem  h  the  subject  of  the  next  section. 

NUMERICAL  ALGORITHM  FOR  THE  MCE  RECONSTRUCTION  OF  THE  EXTINC¬ 
TION  COEFFICIENT 

12.  Although  the  MCE  reconstruction  o f  a*(r)  kaa  beam  conceptually  formulated  aa  a  continuous 
nonlinear  programming  problem,  it  would  ba  emphasised  tkat  from  a  partly  computational  poiat  of 
view,  it  k  aoceasary  to  diacratiaa  tka  problem  aad  rad  oca  it  to  tka  aolutioa  of  a  aoaliaaar  proframmiag 
problem  daftaad  over  a  finite-dimensional  apaca.  Tint,  tkk  raqmiree  a  fiaita-dimenaioaal  approximation 
for  tka  extiactioa  fraction  (wkick  k  aa  lalaita-dimaariraal  paraaaatar).  lb  tkk  aad,  tka  fiaal  eatimata 
q(r)  ia  approximated  bp  a  piacawiaa  cooat  aat  faactioa  of  tka  form 

» 

<(f)  - 

tal 

where  ft  a  •  D  -»  {0, 1}  denotes  tka  ckaractarktk  (iadicator)  faactioa  of  tka  act  S  ia  tka  domain  D 
(rapport  of  extiactioa  fraction).  /Iota  tkat  fta(r)  kaa  valas  1  if  i  £  5  aad  raise  0  otkarwiaa.  It  should 
ba  emphasised  tkat  m  ?n  tka  pracadiag  expreaaim  caa  ba  interpreted  a a  aa  eatimata  for  tka  average 
of  tka  extiactioa  fraction  over  tka  i-tk  range  cell  wkick  k  centered  at  (r<_i  +  u) jl.  Similarity,  a 
piscewiss-conataat  approximation  k  adopted  for  tka  initial  eatimata  p{r),  La. 

#r 

f(r)  *  J2  .r.»(r)- 

Ia  writing  the  pracadiag  expresrions,  it  k  tacitly  aaaamad  tkat  u  >  r<_j  fori*  aad  tkat  tka 

eeqnamcr  of  radial  knots  {u)fml  coaatitataa  a  complete  partition  of  tka  domain  D  (La.,  of  tka  rapport 
of  the  extiactioa  fraction). 

lb.  With  tkeee  approximation*,  tka  discrete  aaalog  of  tka  croaa-entropy  functional  caa  ba  written 

aa  1 

N 

where  if  =  (qu  qq, q*)T  aad  p  =  {puPt,—,Pn)T  are  rector*  ia  tka  Euclidean  apaca  of  dimension 
N  whose  component*  define  tka  piecawi*a-con*taat  approximations  for  q{r)  aad  p(r),  reepectively.  In 
addition,  A r<  a  u  -  r<_i  {»  *  1, 2, . . . ,  K)  denotes  tka  «-tk  incremental  interval  of  tka  range  r.  In  view 
of  tak,  the  dkcretised  form  of  tka  lidar  equation  (cf.  eq*.  (la)  aad  (lc))  reads  aa 

ft  m  5(r<)  -  Ia[*(v<)]  +  * la[«]  -  2 £  Ary.  (4) 

Furthermore,  tka  faactioa  x(r)  aa  well  aa  k  k  implicitly  assumed  to  ba  known  from  tka  functional 
relationship  between  tka  bachsc uttering  aad  extinction  coefficient*.  Now,  combining  *11  tkk,  leads  to 
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tha  wfenwlitiw  of  the  functional  aoaliasar  mraimisatlna  problem  of  eqs.  (3a)  ud  (3b)  as  the  following 
finite  <B—Ail  aoaliasar  minimisation  problem  for  the  determination  of  Um  components  of  the  v sc  tor 

jr 

■?*  IZ*1*  [*/»**] Ar<  M 

*  <«t 

snbjeei  to  Um  constraint 


whart  jfy  is  tha  j-th  observed  lidar  datum  (cf.  aq.  (lb))  wad  jy  is  tha  y-th  medal  (thaontkal)  lidar 
datam  cc-iputed  as  par  aq.  (4). 

14.  It  is  coavsniaat  at  this  poiat  to  mahs  a  few  remarks  concerning  tha  choica  of  tha  number 
of  disenta  intervals  (N)  to  asa  ia  tha  piacawiaa  constant  approximation  of  q{r).  Siaca  tha  optimum 
ran  fa  resolution  ia  tha  direct  ion  of  tha  lidar  baam  path  it  dstarmiatd  by  tha  paisa  duration  r  of  the 
laser  source,  it  is  sensible  to  associate  tha  discretisation  of  q(r)  with  this  quantity.  Consequently,  since 
tha  resolution  call  (or  scatter  call)  has  lenftk  er/2,  it  is  desirable  to  choose  tha  discretisation  so  that 

*  er/2  for  »  «■  1, 2, . . . ,  N.  With  ti  is  schema,  than  now  exists  a  one-to-one  correspondence  between 
tha  estimate  &  at  the  extinction  function  for  tha  »-th  range  call  and  tha  measured  lidar  return  datum  Vi 
from  this  cell.  Hence,  for  this  discretisation  procedure,  the  number  of  discrete  intervals  (N)  chosen  for 
the  pieceerise-constant  representation  of  q(r)  is  equal  to  the  number  of  lidar  data  (M)  employed  in  the 
inversion. 

15.  Although  no  a  priori  constraints  or  bounds  have  been  explicitly  imposed  on  the  values  of  the 
components  of  tha  final  estimate  vector  q,  it  ia  important  to  remark  that,  physically,  these  components 
must  necessarily  be  non-negative.  Strictly  speaking,  tha  non-negativity  constraint  bn  tha  solution  should 
be  imposed  in  addition  to  the  data  constraint  contained  in  eq.  (5b).  However,  the  non-negativity  con¬ 
straint  does  not  need  to  be  explicitly  incorporated  into  tha  formulation,  since  it  is  necessarily  taken 
care  of  automatically  by  virtue  of  the  presence  of  the  logarithm  in  the  cross-entropy  functional.  Now,  it 
should  be  noted  that  eqs.  (5a)  and  (5b)  poses  a  convex  nonlinear  programming  problem  which  can  be 
solved  using  the  ellipsoid  algorithm.  The  ellipsoid  algorithm,  which  gained  prominence  when  Khachian 
[16]  utilised  it  to  demonstrate  the  existence  of  a  polynomial-time  algorithm  for  the  solution  of  linear 
programming  problems,  constitutes  a  simple  and  extremely  robust  procedure  for  solving  convex  program¬ 
ming  problems  such  as  that  presented  by  eqs.  (5a)  and  (5b).  Tbi  procedure  has  the  advantage  of  not 
requiring  the  introduction  of  slack  variables  and  Lagrangian  multipliers  that  has,  to  date,  characterised 
the  numerical  solution  of  all  MCE  problems  (12). 

16.  The  theoretical  basis  for  the  ellipsoid  algorithm  ca  weD  as  an  analysis  of  its  convergence  prop¬ 
erties  can  be  found  in  the  papers  by  Shor  [17]  and  by  Shor  and  Genhovich  [18]  as  well  as  in  the  survey 
by  Bland,  Coldfarb  and  Todd  [19].  A  brief  description  of  the  algorithm,  with  particular  emphasis  on 
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the  details  cotMtiiif  ita  application  to  th«  solution  of  the  kCE  moastructioa  probUm  u  embodied 
ia  sqs.  (5a)  aad  (5b)  foikara.  As  tha  mjm  implisa,  tha  ellipsoid  algorithm  begins  with  the  (pacification 
of  aa  iaitial  tllipooid  So  that  i*  known  to  contain  tha  optimal  aolatioa  to  tha  non-linear  programming 
problam  and  thaa  procaada  to  tha  prescription  of  a  mla  for  tha  generation  of  a  sequence  of  luccejuivel' 
amallar  allipaoida  (La.,  ellipsoids  that  posses*  tm  altar  volumes)  each  of  which  ia  gaarantaad  to  contain 
tha  optimal  aolatioa.  Tha  imptaaostatioa  of  tha  atlipaoid  algorithm  for  tha  aolatioa  of  tha  problam  of 
eqe.  (5a)  and  (Sb)  coaaiata  of  tha  folkariag  thraa  at ape 

17.  SUf  0  :  Giron  tha  iaitial  (prior)  aatimata  p  for  tha  extinction  fraction  aad  a  targat  ralaa  A 
for  tha  miafit,  (apply  iaitial  gmaaaaa  for  tha  final  (poatarior)  estimate  of  tha  extinction  and  for  tha 
initial  allipaoid  matrix  B«.  With  thaaa  iapata,  it  ahoald  parhaps  ba  aotad  that  tha  form  of  tha  initial 
allipaoid  Eo  ia  datarmiaad  as 

Eo  -  e  R" :  (f-  f°)  5 1}, 

whara  ia  tha  cantor  of  tha  allipaoid  and  Ro  it  tha  allipaoid  matrix  which  ia  nacaasarily  real,  symmetric 
aad  positive  definite.  Ia  tha  abaeaca  of  any  prior  information  other  thaa  that  embodied  in  tha  initial 
estimate  p,  it  is  convenient  to  chooaa  f4  «p.  If  lower  aad  apper  bosnda  can  ba  sat  os  tha  extinction 
function,  than  Ro  can  ba  choaaa  so  that  tha  associated  ellipsoid  Eo  ia  the  smallest  ellipsoid  which  encloses 
these  bounds.  Sat  tha  iteration  index  k  *  0.  Steps  1  and  2  that  follow  generate  a  sequence  at  points  qk 
aad  matrices  5t*  that  determine  a  sequence  of  ellipsoids  £*  C  Rr*  for  k  m  1,2,.. .. 

18.  St*p  l  :  Compute  tha  miafit  x3  at  tha  point  <f*.  If  x3  violates  tha  constraint  of  aq.  (5b)  (i.e., 
X3  >  Xm)i  sat  3  m  ^fX3(«*)i  yi*.  tha  gradient  of  tha  misfit  with  respect  to  the  parameter  vector  q 
evaluated  at  j*.  Otherwise,  s«t  Vf8{qk,p),  tha  gradient  vector  of  the  otossk 'tropy  with  respect 
to  f  evaluated  at  qk.  If  aach  component  of  el  ia  lass  thaa  soma  prescribed  tolerance,  then  terminate  the 
procedure  with  qk  as  the  optimal  point;  otherwise,  proceed  to  Step  2. 

19.  Sup  t :  Compute 

dm _ JhL= 

provided  the  quantity  in  the  denominator  is  positive.  Otherwise,  terminate  the  procedure  with  qh  as 
the  best  point  found  by  the  algorithm.  The  updates  for  the  ellipsoid  center  qk  and  the  ellipsoid  matrix 

R*  a re  given  by 

f*41  »f*+<r/(w+i) 

and 

Re+i^-^-fn* - —dr], 

N+l  j» 

respectively.  Replace  k  by  k  +  1  and  return  to  Step  1. 

20.  The  flowchart  for  the  ellipsoid  algorithm  a a  it  applies  to  the  MCE  reconstruction  of  aa  extinction 
function  is  shown  in  Fig.  2.  The  algorithm  is  very  simple  and  can  be  readily  implemented  on  a  computer. 
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Aa  fanplea leatatioa,  which  reeidce  m  *  Microcoft  C  program,  ie  provided  ia  Appendix  B.  A»  a  tad  note, 
it  ia  important  to  potat  oat  that  although  ike  coaetratat  of  eq.  (Sh)  reetrkte  tko  notation  to  Bo  oa  or 
within  rao  donad  hyper-volame  ia  R*,  tko  opitaul  poiat  wkkk  attain*  iaea  tko  croea  iitropy  ia  tkii 
region  laaot  neeueearily  Ik  oa  tko  Mffaeo  of  tko  kypar-eotaam  (do.,  oa  x*  -  Xy)-  Tko  detail*  for  tko 
derivation  of  tkia  naait  an  pm  ta  Appendix  A.  Tkk  (act  kaa  aot  booa  explicitly  incorporated  ta  tko 
iaipkoioatatioa  of  tko  otlipaoid  algorithm  and,  con—queatfr,  caa  ka  aaad  aa  aa  iadopoadoat  chock  of  tko 
optimal  aotatioa  provided  by  tko  algorithm. 

SOME  MCE  INVERSION  EXAMPLES 

21.  Tko  noaka  of  aoato  toata  of  tko  MCE  lid  or  tavenioa  procodan  bjr  ateaaa  of  compator  nmaiatioao 
of  aotay  lidar  data  an  pnaoatad  ia  tkia  taction.  Tko  oatiactioa  pvoftiaa  aaod  ta  tko  ttatalatioaa  an 
compoeitee  kaih  op  from  generic  or  kaa,  raral  aad  maritime  modoia  of  aoroaol  cloada  compikd  bp  tko 
Air  Force  Ceophyeaea  Laboratory  (AFGL)  [20).  Tko  model*  aaod  ia  tko  pnooot  otady  an  tko  Min*  aa 
thoee  aaad  by  BUaoaaotto  [2l|  ta  hi*  aaaoaoaaoat  of  Kiett’*  aaalytfeal  tavaraioa  of  tko  *tagl*-*cau*ring 
lidar  oqaatioa.  All  of  tko  aaroaol  cloada  atiliaod  ta  tkia  atady  bad  a  tkkkaooa  of  1  km  aad  for  rack 
preop eeiflod  extinction  proftlo,  aoiao  free  lidar  data  P(r)  won  |  aerated  for  tweaty  value*  of  rang*  r 
(panning  tko  extent  of  tko  cioad.  Bafon  each  of  them  data  aota  won  i averted  by  tko  MCE  procedure, 
they  won  Irat  corrupted  wilk  Gaamiaa  raadom  aoiao  at  a  prescribed  level  ta  order  to  eimulat*  proceaa 
and  exporimoatal  error*.  Tkia  eorrapted  data  eren  then  logarithmically  traaoformed  to  form  the  S(r) 
data  whick  arm  aaad  ae  ike  iapat  for  Ike  tavonioa  procedan. 

22.  TVe  ftrat  etamlatioa  model  coiiiUea  a  geuark  maritime  *  oool  cioad  at  M  percent  relative  he- 
midity  wilk  a  triaafalar  extinction  prodle  *«(r)  ta  iliaatratod  by  Ike  tolid  line  ta  Fig  3.  The  bee kec alter 
coefficient  f(r)  ta  tkia  example  ia  tmaaiil  to  be  related  to  «»(r)  according  to  0{r)  -  0  0Wo.fr).  Aa  tke 
firet  example,  tereaty  aim  elated  Ildar  data  P{r)  wen  generated  f.  r  tkia  model  and  tkeoe  eo  tee-free  data 
wore  logaritkmkaBy  tranaformed  to  eerve  ae  tke  tapet  for  tke  Kiett  and  tke  MCE  iavereion  procedure 
Klett'e  exact  iavereion  procedure,  with  tke  far-end  boundary  condition  ckooen  to  be  tke  tree  a,  value 
of  1.0  km“‘,  yielded  a  reeelt  tkat  coincided  exactly  eritk  tke  tree  extinction  function  (cf.  Fig.  3).  The 
MCE  tavereion  eolation  with  tke  prior  aatimaie  p(r)  of  ike  extiactioa  function  ckooen  to  bo  a  conetant 
(uniform)  function  oa  tko  rapport  eat  D  (i.e.,  p(r)  -  3.*p©(r))  aloo  gave  a  reeelt  tkat  exactly  matched 

the  tree  extinction  profile.  Fsr  tkia  aoiao- free  data  raee,  it  ia  aecemary  to  eet  <r?  «  1.0  for  •  »  1,2 . Af 

and  to  ehooee  tke  target  miefit  value  Xy  *a  0.001.  Although  both  tke  Kiett  end  tke  MCE  inveteioa  pro¬ 
cedure*  returned  tke  true  extinction  function  under  totae-fre*  condition*,  it  tkoeld  be  empkaeiaed  that 
•  he  MCE  eolation,  anlik*  tke  Kiett  eolation,  did  aot  require  tke  epeciCcation  of  the  correct  boundary 
condition  at  the  far-end  of  tke  aeroeol  clo^d.  Ia  practical  applicatione,  tkia  farmed  boundary  condition 
cannot  be  determined  from  tke  (ivea  data  aad  would  require  aome  additional  infomation  concerning  the 
nature  of  the  aeroeol  cloud.  However,  R  ia  known  [3,21j  tkat  tke  accuracy  of  extiactioa  profile  cecov- 
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end  hy  Um  exact  far-ead  lidar  iiwn in  proeadan  depeade  critically  oa  the  paper  pnacriptin  of  thia 
bewadary  coaditio*.  A a  aa  example,  the  dotted-daahed  Um  ta  Fig.  t  thaws  a  laeoaatractad  extiactioa 
faactioa  (for  the  maim  free  dau  oat)  mU|  KleU'a  proeadan  with  the  farmed  handary  caaditioa  apari- 
ftad  aa  a,  •  10  km-1.  Ohaarva  that  althoayk  the  aolatioa  coavnyae  ta  the  tnt  axtiactioa  coeffkimt  at 
the  naread  of  the  cload,  than  it,  aanrthalsea,  a  aaoderate  error  hi  the  datanaiaatioa  of  the  axtiactioa 
at  the  far  aad  of  the  dead.  Aa  a  faal  poiat,  it  ahoatd  ba  aotad  that  akhoaghths  MCE  proeadan  doaa 
aot  raqaira  tha  apadicatioa  of  aa  appropriate  far-ead  boaadary  niae,  it  haa  da  diaadnatay*  of  bain* 
con  petal  ioaally  ton  expo  tain  that  tha  Klett  procadara. 

33.  Tha  text  example  coaddara  tha  ten  raaliatic  cam  of  toiay  Hdar  date.  Flf.  4  compana  tha 
raatlt  of  tha  MCE  aolatioa  with  Klatt'a  aaalytk  aolatioa  for  tha  cam  of  bdar  data coataadaaled  with  l.S 
parent  RMS  Cataaiat  toiaa.  For  thia  example,  tha  MCE  aolatioa  waa  coaapatad  with  tha  prior  eatimate 
p(r)  chooaa  to  ba  a  coaotaat  ftactioa  aad  tha  tldi  taryat  vdee  xV«aa  talO.  Note  that  tha  MCE 
aolatioa  ie  eaaooth  aad  approximates  tha  tnaa  axtiactioa  profit  nry  waft  Oa  the  other  haad,  tha  Kletl 
aolatioa  which  waa  obtaiaad  with  tha  correct  epecilcaiioa  of  the  f tread  bond  cry  vdee,  diaplayad  a 
moderate  amoeat  of  oocillatioa  at  tha  far-ead  of  tha  cload.  Obearn  that  thieeadllatioa  it  aot  praaaat 
at  tha  neer-ead  of  the  cload  aad,  ia  thia  reyioa,  the  approxhaatioa  of  tha  Khtt  aolatioa  to  the  tree 
extiactioa  profit  ia  aa  yood  aa  that  provided  by  the  MCE  aolatioa.  For  am  a  email  lent  of  aoiaa, 
the  MCE  ianraioa  procadara  provides  aa  iatriaaicaUy  aaeoothar  aolatioa  thaa  that  gima  by  the  Klett 
iavenioa  procadara.  AHhoagh  both  ianraioa  procedaraa  an  at  able  with  nee  pact  to  partarbeiioM  ia  the 
data,  the  MCE  ianraioa  ia  aot  aa  greatly  affected  by  the  aoiaa  aa  it  the  Klett  ianraioa. 

34.  Next,  eome  axampiee  which  illaatrate  certaia  propertiaa  of  tha  MCE  iaeanioa  proeadan  with 
reap  act  to  the  choke  of  the  prior  eat  ha  ate,  the  choice  of  tha  atiaft  taryat  taka,  aad  the  lent  of  aoiaa 
comptiag  tha  Udar  data  an  pnoeated.  Tha  reealt  of  aa  MCE  ianraioa  for  •  t  parent  lent  of  RMS 
Gaaaaiaa  aoiaa  ia  ahowa  ia  Fly.  S.  Tha  dotted  Um  ia  Fig.  I  cpmepoadt  to  aa  MCE  ianraioa  with  tha 
prior  eatimate  p(r)  ckoen  to  ba  a  coaataat  or  aaifona  faactioa.  Agaia,  tha  atiaft  taryat  valae  xlamto  WM 
eat  to  30.  Siace  tha  prior  eatimate  ia  eatfona,  tha  MCE  iannioa  reealt  ia  (haply  tha  maxima  ia  mtropy 
raatoratioa  of  the  axtiactioa  profile  aad  aa  rock,  provides  the  moet  coassrvetiee  (eaiform)  eatimate  of 
tha  axtiactioa  coasisteat  with  tha  fine  data.  Obaam  that  tha  iavardoa  waa  take  aacceadal  eiaea  tbe 
rank  y(rl  (i.a.,  the  fad  or  poeterior  eatimate  of  the  extiactioa  faactioa)  ckmty  approximates  tha  true 
ntiactioa  profile  a,(r),  evea  ia  tha  praaaaca  of  dmeietad  aoiaa.  Next,  tha  target  mkifit  valae  waa  eat 
to  33  with  ao  ekaagee  la  the  other  parameters.  The  MCE  ianraioa  reealt  far  thia  case  ia  displayed  by 
tha  dashed  Uaa  ia  Fif .  $.  Slate  the  aoiee  level  wae  slightly  overestimated  la  tkia  example,  obearn  that 
the  raealtiay  ianraioa  ia  ovenmoothed  (via.,  the  peak  of  tha  extiactioa  faactiax  ia  aadenetimsted)  aad, 
iadaad,  the  fad  estimate  ia  thia  cam  more  cloeely  reeembiea  the  i  ait  id  eaifanB  eatimate  thaa  doee  the 
fiad  estimate  ia  the  previoxa  example.  Tbit  ehoeid  aot  be  too  earpriaixg  dace  ia  the  cam  Xu  —  oo, 
the  data  does  aot  coastraia  the  aolatioa  at  dl  and  tha  miaimam  cross  a atropy  aolatioa  is  achieved  whea 
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*(r)  m  p(r),  vi>.  when  the  laal  dim  ate  to  itotkil  to  Dm  initial  eetimate.  Thto  to  perfectly  utml 
timet  one  woald  M  Mr*  tin  tail  din  ate  to  difer  froan  the  initial  dimate  *tot  Um  ototmd  data 
coaatota  of  pwe  mom.  b  tooald  bo  murtod  ttoi  to  achieve  a  food  inversion,  to  to  aecamary  to  oxarctoo 
jadgoaeeai  ia  tto  choke  of  too  audit  larfot  vahte,  via.  toto  nht  ahoald  bo  eboaoa  ao  toot  too  eotatioa 
Ito  too  gtoaa  data  within  aeaao  exported  or  prescribed  toiaraaca  that  edeqaaiety  refeete  too  severity  of 
too  aotoa  tempting  too  data. 

38.  Tho  aocoad  aiaanlaiioa  aaodai  considered  bar*  to  coastrocted  from  a  coaapooito  coafigaratioa 
of  too  generic  arbaa,  roraJ  aad  maritime  atroooto  at  70  relative  bom  id  tty  with  a  raaabiag  extiactioa 
proftlo  displayed  by  tho  ootid  lino  ia  Fig.  8.  For  tkio  example,  too  faaetioaal  depeadeace  between  tho 
beckec  at  taring  aad  oxtiactioa  coefficient*  mirror  the  regional  aaroooto  oood  ia  the  coiatractioa  of  the 
eoatpooite  model.  Boaca,  the  faaetioaal  depeadeace  to  doocribod  by  a  ptoce-coaataat  reiationahip  of  the 
form 

fi[r)  m  |o.OOMato.o  s*)(r)  +  O.O306ftjo.m,e.ee)(r)  +  0.021p|4.wl.Ooj(r)  jo«(r). 

Twenty  eimaiatod  bdar  data  were  generated  for  toto  exam  pie  aad  eorraptod  with  Qaaeeiaa  aotoe  at  the  8 
percent  RMS  aotoe  keel.  With  the  choice  of  p(r)  ao  a  coaotaat  or  aaiform  fa  action  (i.e.,  p{r)  -  3.8|to(r)) 
and  with  tho  mtofit  target  ealae  eat  at  Xymtt a  80,  the  raealtiag  MCE  recoaatractioa  of  the  oxtiactioa 
profle  to  iOaetratod  by  the  dotted  line  of  Fig.  8.  The  ramltieg  iaetrtioa  to  rtaaoaabie  for  the  elated 
aotoe  level  Tho  daahrd  line  of  Fig.  8  ehowa  the  reeah  of  MCE  iavereion,  bat  thto  time  with  aa  initial 
eetimate  p(r)  given  by 

p(r)  -  (1.0+  I0.0r)^e  t)(r)  +  (11.0-  10.0r)^.itl  ot(r).  (8) 

Obeerve  that  tkio  iavereion  reoab  to  qnlte  aimilar  to  tho  preceding  one  with  the  exception  that  the  peak 
ia  the  extinction  profle  at  r  -  0.8  km  to  (lightly  better  defaed.  Thto  tooald  not  be  too  twrpriiing  line* 
the  prior  eetimate  p{r)  for  thto  cam  explicitly  encoded  the  preeeaco  of  tkio  peak  aad  the  MCE  iavenion 
attempt*  to  chooe*  the  laal  eetimate  to  be  ao  timilar  to  the  iaitial  eetimate  aa  poeeible  while  eatiefying 
the  conet  relate  impoood  by  the  data. 

38.  To  provide  aa  idea  on  the  biea  aad  variability  ia  the  MCE  rec  on  (traction  of  the  extinction 
feactioa,  the  MCE  iavereion  precede re  wee  applied  fifty  timea  to  the  eame  iim elated  (aotoeleea)  data  art 
(i*neritrd  for  the  compoeite  aeroeol  model  deecribed  above)  to  which  fifty  eeta  of  Gaeeeiaa  aotoe  have 
been  added  at  the  8  parent  RMS  level.  For  thto  Monte  Carlo  aimniatioa,  the  biaa  B(r)  aad  variability 
V  (r)  ia  the  MCE  invert  ioe  retell  to  eetimated  aa  foil  owe  : 
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*hn  f'(r)  took*  Sk  MCE  wtowtnrtici  of  the  extinction  function  far  i-tfc  dote  aol  tad 

i(r)“s5T^w 

tel 

io  the  MU  of  the  50  MCE  iavenioue.  The  bias  is  (imply  tke  mu  of  the  fifty  differences  between 
Um  reconstructed  aad  Um  true  extinction  functions.  The  bmu  f(r)  of  the  SO  inversions  io  ikowt  by 
UM  dotted  curve  io  Fig.  7.  Tho  error  bon  (boor  om  standard  deviations  Emit  (U.,  Ik  variability 
V  (r))  computed  from  tbo  SO  hum.  Tbooo  results  indkate  tbot  tbo  mvewiau  proeodon  porforau  quite 
adequately  especially  when  out  roc  (Ik  tbot  tbo  doU  oot  couoioto  ofouly  twenty  ooioy  lidar  measurements. 
Tbo  biao  ia  tbo  recovered  oxtiactiou  profik  io  oa  tbo  order  of  0.2S  km-1  whereas  tbo  varisbifity  io  ou  tbo 
order  of  0.30  to  0.S0  km-1.  Observe  tbot  tbo  variability  ia  tba  MCE  neoaatractioa  of  tbo  oxtiactiou 
coefficient  k  larger  at  tbo  far  ead  of  tbo  cloud  tbaa  at  tba  boot  aad.  Tkia  k  due  to  tbo  fact  that  tbo 
lidar  retaro  aigaal  P[r)  dacroaom  with  raage  dot  to  attenuation;  hence,  tbo  comopoadiag  logarithmic 
raago-adjaatod  aigaal  5(r)  oxbibito  aa  iaeroaaiag  ooko  level  with  increasing  raage. 

27.  Ia  actual  practice,  out  aever  knows  tbo  procka  relatioaakip  bet  woe  a  tbo  backacatteriag  and 
extinction  coefficient*.  To  aunalate  tbooo  circumstances,  tbo  MCE  ioveraka  for  tbo  example  coaaidered 
ia  Fig.  0  k  recomputed,  bat  tbk  tune  with  a  mio-e pacification  ia  tbo  functional  relatiouabip  between  the 
backacatteriag  aad  oxtiactiou  coefficient*,  vk. 

fi(r)  m  |o.01/<fol0.3*|(r)  -f0.02Stqo.3n,i.  Q|W}«eW- 

Tbo  initial  ootimato  p(r)  wa*  choaon  ao  pir  oq.  (0).  In  order  to  partially  account  far  tbo  fact  that  tba 
backacatteriag- to-extinction  ratio  k  poorly  known,  the  misfit  ta*«at  value  xL-so  wa*  to  27.  The 
invenioa  result  k  shown  by  the  dotted  carv*  ia  Fig.  8.  Tbo  similarity  between  the  MCE  recovered  and 
tba  true  extinction  profik  k  aa  good  ao  can  be  expected  especially  ia  light  of  tbo  incorrect  specification 
of  the  backac alter- to-extinction  ratio  aad  the  level  of  noise  ia  the  data.  Observe  that  the  main  peak  in 
tbo  extinction  function  at  the  midpoint  of  the  clovd  k  well-determined.  Tbo  inverted  reeuit  show*  the 
two  smaller  peaks  also,  but  their  values  us  underestimated. 

CONCLUSIONS 

28.  Aa  inform  at  i<u>theore  tic  approach  based  ou  the  priacipk  of  minimum  crons  entropy  (MCE) 
that  allows  for  the  rational  combination  of  prior  lnfornu>taa  on  Um  extinction  profile  with  additional 
information  coatakod  ia  the  aoky  lid  v  return  data,  baa  been  utilked  to  formulate  aa  inversion  procedure 
for  the  riagk-eeatteriag  lidar  eqaation.  Since  noa-saiquriwo  in  the  inversion  of  the  Ildar  equation  ia 
well-known,  tke  MCE  invenioa  procedure  enables  Um  objective  reconstructs.^  of  a  preferred  model  of 
tke  aerosol  cloud  that  not  only  minimises  extraaeons  features  in  tke  solution  not  reqi  .vd  ia  order  to 
interpret  tke  data,  bu«  also  optimally  accounts  for  errors  in  tbs  data.  Indeed,  the  mature  of  the  model 
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Nconnd  depends  as  nmck  cm  tW  diU  wen  Moatk  data  themselves.  A  good  inversion  rmk,  which 
dose  not  nunlead  the  user,  cam  only  b«  obtained  it  ail  ntilabh  iaformaiioa  cm  th<  problem  is  effectively 
■tilined  aad  this,  of  wins,  inchedes  the  information  on  the  level  of  aoiee  corrupting  the  data. 

29.  A  Maple  aad  robnst  aaaxrical  pncedaie  baaed  oa  the  eilipeoid  algorithm  hae  beea  developed 
far  the  eelatioa  of  the  highly  aoaBaeer,  ronetriiaed  auaimisatioa  problem  repaired  to  perform  the  MCE 
recoaetraotioa.  Ifamerkal  raeaha  ohta'nsd  esiag  ehwalated  aoiey  lidar  data  show  that  the  algorithm  for 
MCE  recoastrectioa  provider  stable  aad  well- behaved  eetimatea  for  the  extinction  function.  Ia  addition, 
thia  ammvrical  procedare  seesaa  to  be  qaite  iaaeasitive  to  the  ^-ivea  iaitial  gaeee  for  the  final  estimate 
g(r).  Except  far  the  iaitial  pheaa,  the  tllipeoid  algorithm  geaerally  exhibited  coavergeace  to  the  optimal 
MCE  eolation  at  a  liaear  rate,  a  reeah  that  haa  beea  predicted  by  theory  [17].  Conaeqaeatly,  while  the 
algorithm  displays  quite  a  rapid  coavsTgeace  to  the  aear-optimal  eoletioa  after  the  first  twenty  to  forty 
iterations,  it  nevertheless  exhibits  a  sices  coavergeace  rate  a  ear  the  optimal  solution.  Since  is  take*  from 
aboat  30  to  60  seconds  to  compnte  aa  MCE  inversion  solntioa  oa  aa  IBM  PC-XT,  the  algorithm  in  its 
present  form  is  not  suitable  for  a  real-time  recovery  of  extinction  profiles,  aaleee  either  aa  algorithm  with 
a  better  coavergeace  rate  aear  the  optimal  eolation  is  developed  or  a  computer  with  a  faster  processor 
is  used. 
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Schematic  illustration  of  the  geometric  interpretation  of  the  principle  of  minimum 
cross-entropy.  Here,  the  final  estimate  q  is  the  information-  or  1-projection  of  the 
initial  estimate  p  onto  the  admissible  set  IT  Consequently,  q  is  as  close  to  p  as 
possible  in  the  information  measure  sense  while  at  the  same  time  satisfying  the 

new  information  encoded  in  the  data. 
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Figure  2 

Flowchart  of  the  eltipeoid  algorithm  used  to  eoive  tha  convex  optimization 
problem  required  to  obtain  the  minimum  croea- entropy  recovery  of  the 

extinction  function. 
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Figure  3 

A  comparison  of  the  aerosol  extinction  functions  recovered  by  the  MCE  and  Klett 
inversion  procedures  for  the  caise  of  noise-free  iidar  data.  The  true  extinction 
function  (solid  line)  is  representative  of  a  generic  maritime  aerosol  cloud  at  99% 
relative  humidity.  The  MCE  solution  and  the  Klett  solution,  with  a  correctly 
specified  far-end  boundary  condition,  returned  extinction  profiles  that  coincided 
with  the  true  extinction  coefficient.  The  dashed  line  shows  the  Klett  solution 
with  an  incorrectly  specified  far-end  boundary  value  on  the  aerosol  cloud. 


UNCLASSIFIED 


EXTINCTION  COEFFICIENT  (1/km) 


UNCLASSIFIED 


SM  1225 


Figure  4 

A  comparison  of  tha  aerosol  extinction  functions  recovered  by  the  MCE  and  Klett 
inversion  procedures  for  the  case  of  Ildar  data  contaminated  with  1.5%  RMS 
Gaussian  noise.  The  Klett  solution  was  obtained  with  the  correctly  specified 
far  end  boundary  value  for  the  aerosol  cloud.  Note  that  tha  MCE  solution  is 

much  smoother  than  the  Klett  solution. 
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Figure  5 

An  example  of  minimum  cross-entropy  recovered  aerosol  extinction  functions. 
The  dotted  end  dashed  lines  indicate  MCE  inversions  with  target  misfit  values  of 
20  and  32.  respectively.  The  input  for  these  inversions  consisted  of  20  Ildar  data 
degraded  with  3%  RMS  Gaussian  noise. 
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Figure  6 

Another  example  of  minimum  cross -entropy  recovered  aerosol  extinction  functions. 
The  true  extinction  function  (solid  line)  is  a  composite  constructed  by  piecing 
together  generic  aerosol  cloud  characteristics  from  urban,  rural  and  maritime 
environments.  The  dotted  line  indicates  the  inverted  extinction  function  for  a 
target  misfit  value  of  20  and  a  uniform  initial  estimate.  The  dashed  fine  shows 
the  inverted  extinction  function  using  a  target  misfit  value  of  20  and  the  initial 
estimate  as  the  triangular  extinction  profile  shown  by  the  solid  line  in  Figure  5. 
These  inversions  were  obtained  from  twenty  Udar  data  which  have  been  corrupted 

with  5%  RMS  Gaussian  noise. 
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Figure  7 

Same  example  as  in  Figure  6  but  the  dotted  line  shows  the  means  of  50  MCE 
inversions.  The  error  bars  indicate  the  standard  deviations  /o'*  the  inversions 
computed  from  the  fifty  noise-corrupted  iidar  data  sets.  Each  data  set  consisted 
of  20  data  that  has  been  corrupted  with  5%  RMS  Gaussian  noise. 
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Figure  • 

Same  example  a*  in  Figure  8  but  the  dotted  line  show*  the  recovered  extinction 
profile  for  a  mltspecified  backscatterlngto-extinctlon  ratio.  The  target  misfit 

value  we*  set  to  27  for  the  inversion. 
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APPENDIX  A 

90.  The  MCE  recoastractioa  of  tW  extinction  coeflkieat  of  aerosol  cMi  from  iacompUu  ud 
uurtiii  lidar  rot  on  data  is  formalatad  aa  the  problem  of  fiadiag  tka  particular  aolatioa  tkat  poaaaaaaa 
tk«  smaliset  cross  satropy  relative  to  aoaao  initial  (prior)  nodal  and  tkat  at  tka  mm  timo  provide*  a 
satisfactory  snatch  between  tka  modal  (tkooratkal)  and  actoal  data,  witk  tkk  matck  baiag  qaaatified 
by  a  data  miait  of  tka  form 


It  will  bo  akowa  tkat  tka  aolatioa  of  miaimam  croaa  latropy  liaa  oa  tka  boaadary  of  tka  allowabl*  misfit 
region  dafinad  abort  (vis.,  oa  tka  aarfaca  x*  “  Xu)  ratkar  tkaa  ia  ita  interior.  la  otkar  words,  tka 
aolatioa  of  miaimam  croaa  eatropy  ia  aocomarily  tka  oaa  tkat  corraapoada  to  tka  Urfoat  parmiaaibla 
misfit. 

91.  Tkia  raaalt  is  aaaily  damoaatratad  by  eontrw^-Uon.  Sapposa  tkat  tka  atiabaan  cross  aatropy 
aolatioa  f*(rj  raaalta  ia  tka  data  misfit  x*1  <  Xu,  i.a.  tka  aolatioa  liaa  ia  tka  iatarior  of  tka  misfit 
rtfioa.  Now  coaaidar  a  aolatioa  giraa  by  f(r)  -  *g*(r)  (r  6  (0, 1))  witk  tka  cortespoadiag  data  misfit 
X*.  Next,  express  g*  ia  terms  of  x‘*  by  avbstitatiag  f (y)  iato  #q.  (A.l)  sad  performing  aosaa  simple 
algebraic  maaipalatioas  to  obtaia 

X*  >  x#*  +  (l  -  M)«,(g*(r))  +  (l  -  **)«*(♦»)  +  (laM)«*(»*(r))  +fc*la*(x).  (A2) 

where  C|,  <j,  aad  e*  are  coastaaU  tkat  dapaad  oa  g*(r).  Now,  obsanrs  from  *q.  (A.2)  that  if  x  is 
selected  to  be  sa  flkiently  close  to  oaa,  it  is  possible  for  X*  to  be  strictly  Isas  tkaa  Xl,  slice  x*‘  £  xl 
by  aoppoeitioa  (vis.,  tka  poiat  f(r)  liaa  strictly  ia  tka  iatarior  of  tka  acceptable  misfit  region).  However, 
note  that  tka  valve  of  tka  cross-* atropy  fvactioaai  at  f(r)  is  lass  tkaa  tka  aaaaiaed  miaimam  valaa  of 
the  fvactioaai  at  f*(r)  siaca 

*  (f.  *)  m  *  fo  f*  (f)  1*  [f *  M/fW]  dr  +  v  la  v  g 1 *  (r)  dr 

-*ff(f*,p)  +  *lax/  f*(r)dr. 

Jd 

Consequently,  ia  coatradictioa  witk  tka  original  kypotkaais,  tka  MCE  aolatioa  mast  lit  oa  tka  boaadary 
at  the  allowabl*  misfit  fvfioa. 
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APPENDIX  B 

11  Tko  program  gma  un  tkk  Appcadix  it  u  iapknoUtio*  of  tk  ellipsoid  algorithm  for  tko  MCE 
roeoaotractioa  of  oxtiactio*  futtiou  for  mtooo!  dork  It  k  writtoa  ia  Microsoft  C  which  coaforms 
to  tko  lorthcomiag  America*  Natioaal  Staadard*  Iaotitato  (ANSI)  standard  for  tko  C  laagaago.  It 
skoald  bo  mpkmixd  tkat  tko  program  km  boo*  writtoa  largely  for  clarity  ratkor  tkaa  (pood  aad,  as 
rack,  skoald  aot  bo  ragardod  as  a  prodoctioa  softwaro  impkmoatatioa  of  tko  algaritkm.  It  k  primarily 
iatoadod  for  rosoarck  parpooso. 
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# include  <stdio.h> 
# include  <math.h> 
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/* 

/* 

/* 

/* 


***  / 


_ V 

- V 

V 

Title  :  An  Information-Theoretic  Method  for  the  Inversion  of  */ 

the  Lidar  Equation  */ 

V 

Date:  October/  1987  */ 

■  V 

Author:  Eugene  Yee  */ 

Defence  Research  Establishment  Suf field  */ 

Ralston/  Alberta  TOJ  2N0  */ 

*/ 

- */ 

- - - - v 

*/ 

Name  -  ellipsoid(N,N  DATA,  eps,chi_sqrm,max_iter,  data,  */ 

variance, R_pos, vie  r,delta_r,x,y,rmatrix)  */ 

*/ 

Purpose  -  To  apply  the  ellipsoid  algorithm  to  solve  the  */ 

nonlinear  programming  problem  required  to  compute  */ 

the  minimum  cross-entropy  reconstruction  of  the  */ 

extinction  coefficient  functions.  */ 

V 

- v 

Description  of  parameters  in  the  argument  list:  V 

V 

N  number  of  variables  */ 

MDATA  number  of  data  points  */ 

eps  accuracy  required  in  each  element  of  gradient  */ 

vector  V 

chisqrm  target  value  for  the  misfit  */ 

maxitor  maximum  allowable  number  of  iterations  of  */ 

ellipsoid  algorithm  */ 

data  array  that  stores  the  lidar  return  data  */ 

variance  variance  in  the  values  in  the  data[]  array  */ 

Rpos  array  storing  radial  values  r  at  which  lidar  */ 

data  stored  in  data[]  are  measured  */ 

vec  r  array  storing  radial  values  r  at  which  */ 

extinction  function  will  be  estimated  *>/ 

deltar  array  that  stores  incremental  radial  values  */ 

between  elements  specified  in  array  vec  r[]  */ 

x  stores  current  approximation  of  extinction  */ 

function.  An  initial  approximation  should  be  */ 

provided  on  entry,  which  will  be  replaced  by  */ 

the  MCE  estimate  on  exit.  */ 

y  stores  the  initial  (prior)  estimate  of  the  */ 

extinction  function  */ 

rmatrix  used  to  store  the  currant  ellipsoid  matrix  */ 

generated.  An  initial  ellipsoid  matrix  should  */ 

be  provided  on  entry.  */ 

*/ 

- - - - - - */ 

Remarks: 


1)  The  user  must  provide  two  functions  : 


/*  (i)  double  model_data (N ,  Rj , veer , del tar ,  x)  */ 

/*  which,  given  the  range  Rj  and  the  parameter  vector  x[]  V 

/*  calculates  the  corresponding  model  or  theoretical  V 

/*  lidar  datum.  */ 

/*  (ii)  grad_dj (N,Rj ,vec_r,delta_r,x,gnd_dj)  V 

/*  which,  given  the  range  Rj  and  x,  calculates  the  gradient  */ 

/*  of  the  associated  model  datum  with  respect  to  the  para-  */ 

/*  meter  vector  x  (current  estimate  for  extinction  function)  */ 

/*  and  places  the  result  in  the  array  gnd_dj[].  V 

/*  */ 

/*  2)  The  ellipsoid  matrix  rmatrix()n  is  assumed  to  be  */ 

/*  dimensioned  in  the  calling  program  with  column  */ 

/*  dimension  50.  V 

/*  V 

/*  _  V 

/* - */ 

/a******************************************************************/ 


ell ipsoid (N , NDATA , eps , chi_sqrm , mnx_iter , data , variance , Rpos , veer , 
del ta_r , x , y , rmatr ix ) 

unsigned  int  N,N_DATA,max_iter; 
double  eps,chi_sqrm; 

double  dataf ] , variance [ ] ,R_pos[ ] ,vec_r( ) ,delta_r[ J ,x() ,y [ ) ; 
double  rmatrix(] [50] ;  ~ 

( 

register  int  anlnt,anlntl; 
unsigned  int  iter_count; 

double  chi_sqr , tmp , tmpl , gl , Rj , FACTOR , NP_ONE  t 
double  gnddj [50] ,grad[50] , tmparr[50] ; 
double  model_data ( ) ; 
double  euclid_norm() ; 


iter_count  «  0; 

FACTOR  -  N*N/(N*N  -  1.0); 
HP  ONE  -  N  +  1.0; 


/*  »»»»»»»»»»>»»»»»»»»>»»»»>  */ 
/*  Begin  by  computing  value  of  constraint  function  */ 
/*  (i.e.,  the  misfit)  at  the  currant  value  of  the  */ 

/*  parameter  vector  x.  */ 

/*  But  first,  test  to  see  that  maximum  number  of  */ 
/*  specified  iterations  of  the  algorithm  has  not  */ 
/*  been  exceeded.  */ 

/*  »>»»»»»»»»»»»>»»>»»>»»»>»»  */ 

iter_count  +■  1;  /*  Count  no.  iterations.  * / 

if (iter_count  >  max_iter) 

( 

printf{"  Maximum  no.  iterations  exceeded. \n") ; 

return ; 


) 


/*  »»»»»»»»»»»»»»»»»»»»»»»  */ 
/*  Now  begin  the  computational  work  of  STEP  1.  */ 

/*  »»»»»»»>»»»»»»»»»»»»»»»>  */ 


chi_aqr  -  0.0; 

for(anInt  -  0;  anlnt  <  N_DATA;  anlnt++) 

{ 

Rj«R_poe[ anlnt] ; 

tmp*data[ anlnt ]  -!-model_data  (N ,  R j ,  vec_r ,  delta_r , x) ; 
tmpl»tmp*tmp; 

chi_sqr  +■*  tmpl/variance  [ anlnt  ] ; 

) 

gl  -  chi_sqr  -  chi_sqnn; 

if(gl  >  eps)  /*  Is  constraint  violated  ?  */ 

{ 


/*  »»»»>»»»»»»»»»»»»»»>»»»».  */ 
/*  Calculate  the  gradient  of  the  constraint  V 

/*  function.  */ 

/*  »»»»»»»»»»»»»»»»»»»»»»»  */ 


for(anInt  -  0;  anlnt  <  N;  anlnt++) 

grad ( anlnt ] *0. 0;  /*  Zero  gradient  array.  */ 

for(anInt  -  0;  anlnt  <  N  DATA;  anlnt++) 

( 

Rj»R_P°s [ anlnt ] ; 

graddj (N,Rj ,vec_r,delta_r,x,gnd  dj); 

tmp  -  data[anInt]-model_data(N,R5,vec_r,delta_r,x)  ; 

tmpl  »  tmp/variance ( anlnt ] ; 

for (anlnt 1  ■  0;  anlnt 1  <  N;  anlnt 1++) 

grad[anlntl]  —  2.0*tmpi-*gnd_dj  [anlntl] ; 

) 

) 

else 

< 


/* 

/* 

/*, 

/* 


Calculate  the  gradient  of  the  objective 
function  (i.e.,  the  cross-entropy  ^unction) 
»»»»>»»»»»»»»»»>»>>»»»»»»> 


V 

V 

V 
*/ 


for(anInt  -  0;  anlnt  <  N;  anlnt++) 

grad(anlnt]  -  (1.0+log(x[anlnt]/y(anlnt]) ) ; 


/*  »»»»»»»»»»>»>>»»>»»>»»»»»»  */ 
/*•  Test  for  convergence  */ 

/*  >>>>»>»>»>»>>>»>>>>>>>>>>>>>>>>>»>>>>>>»  */ 


for (anlnt  -  0;  anlnt  <  N;  anlnt++) 

if (gradfanlr.t]  >  eps)  goto  steptwo; 
return;  /*  Convergence  of  algorithm:  exit  */ 


/* 


*/ 


/*  >»»»»»»»»»»»»>»»»»»»»»»»»  */ 
/*  Normalize  the  gradient  array  grad[]  by  its  */ 

/*  length.  The  Euclidean  length  of  grad[]  is  */ 

/*  computed  using  the  function  euclid_norm() .  */ 

/*  »»»»»»»»»»»»»»»»»»»»»»»»  */ 


tmp  «  euclid_norm(N,grad) ; 
for(anInt  -  0;  anlnt  <  N;  anlnt+t) 
grad[anlnt]  /-  tmp? 


/*  »»»»»»»»»»»»»»»»»»»»»»»»  */ 

/*  In  this  first  portion  of  STEP  2,  we  make  all  */ 

/*  the  needed  calculations  for  the  update  of  the  */ 

/*  parameter  vector  (which  determines  the  center  */ 

/*  of  the  ellipsoid)  and  the  ellipsoid  matrix  */ 

/*  rmatrix[)[3.  V 

/*  >»»»»»».>»»»»»»»»»»»»»»»»>  */ 

for  (anlnt  ■  o?  anlnt  <  N;  anlnt+-f) 

( 

tmparr [ anlnt ] »0 . 0 ; 

for (anlnt 1  -  0;  anlntl  <  H?  anlntl++) 

tmparr[ anlnt}  +-  rmatrix[anlntnanlntl]*grad(anlntl] 

tmp-0 . 0 ; 

for  (anl.it  -  0;  anlnt  <  N?  anlnt++) 

tmp  +-  grad [anlnt ]*tmparr[ anlnt } ; 
if (tmp  <•  eps) 

( 

printf ("  Ellipsoid  no  longer  positive  definite. \n") ; 
return;  /*  exit  */ 

) 

tmpl  -  sgrt(tmp); 


for (anlnt  -  0;  anlnt  <  N;  anlnt++) 

tmparrfanlnt]  ■  -tmparrfanlnt J/tmpl; 


/* 

/* 

/* 

/* 

/* 


After  these  preliminary  computations,  we  are 
now  ready  to  update  the  paramter  vector  x[] 
and  the  ellipsoid  matrix  rmatrix[][]. 


*/ 

*/ 

*/ 

*/ 

*/ 


tmp  -  2 . 0/NP_ONE  ? 

for (anlnt  ■  0;  anlnt  <  N;  anlnt++) 

{  /*  First  update  parameter  vector.  ,  */ 

x[anlnt]  +»  tmparr[anInt]/NP  ONE? 
for(anlntl  -  0?  anlntl  <  N?  anlntl++) 

{  /*  Next  update  the  ellipsoid  matrix.  */ 

rmatrix[ anlnt] [anlntl]  — 
tmp*tmparr [ anlnt ] *tmparr[ anlntl] ? 
rmatrixf anlnt] [anlntl]  *«  FACTOR; 

) 

) 


/*  »»»»»»»»»»»»>»: 


/*  After  thesa  updates,  go  back  to  STEP  1  and  */ 

/*  start  all  over  again.  V 

/*  »»»»»»»»»»»»»»»»»»»»»>»»  */ 

goto  step_one; 

} 

/♦a*****************************************************************/ 


/*  */ 

/*  Name  -  grad_dj (N,Rj ,vec_r,dalta_r,x,gnd_dj)  */ 

/*  */ 

/*  Purpose  -  To  compute  the  gradient  of  the  j-th  model  data  */ 

/*  point  d[j]  with  respect  to  the  parameter  vector  */ 

/*  x  -  (q[l] , . . . ,q[N] )  with  the  result  placed  in  the  */ 

/*  gnd_dj  [  ] .  Rj  is  the  range  corresponding  to  the  */ 

/*  model  datum.  V 

/*  V 


/ ************************************ ************************* ******/ 
gradd j (N , R j , vec_r , del t a_r , x , gnd_d j ) 
unsigned  int  N; 

double  Rj,vec_r[],delta_r[],x[],gnd_djC]; 

{ 

register  int  anlnt; 
double  K; 

K  -  1.0; 

for (anlnt  -  0;  anlnt  <  N;  anlnt ++) 

( 

if ((anlnt  -»  0)  &&  (vec_r[anlnt]  >Rj)) 

gnddj [anlnt]  »-2.0*delta_r [anlnt]  +  K/x [ anlnt ] ; 
else  if (vec_r( anlnt]  <  Rj) 

gnddj [anlnt]»-2.0*delta_r[anlntj ; 
else  if (vec_r(anlnt-l]  <-  Rj  &&  Rj  <  vec_r[anlnt] ) 

gnd_dj (anlnt ]«-2.0*delta_r[anInt]  +  K/x[anlnt]; 
else  gnd_dj [anlnt ] *0.0; 

) 

) 


/A******************************************************************/ 


/*  V 

/*  Name  -  model_data(N,Rj,vec_r,delta_r,x)  */ 

/*  */ 

/*  Purpose  -  To  compute  the  j-th  component  of  the  model  data  */ 

/*  vector  corresponding  to  range  Rj.  */ 

/*  V 


«  / ************  ************************************************** *****/ 
double  model_data(N,Rj,vec_r,delta_r,x) 
unsigned  int  N; 

double  Rj ,vec_r[] , delta  r[],x[]; 

{ 

register  int  anlnt; 

double  trap , CKST1 , CKST2 , CKST3 , K; 


tmp»0 . 0  ? 

CKST1  -  0.01; 


CKST2  -  0.025; 
CKST3  -  0.025; 
K  -  1.0; 


for  mint 
{ 


0;  anlnt  <  N;  anlnt++) 


if ((anlnt  —  0)  44  (vec_r[anln 
tap  +»  x[anlnt]*delta_r[a 
•lse  if (vec_r [anlnt]  <-  Rj) 

tap  +«  x[anlr.t]*delta_r[a 
•Isa  if ( (vec_r[anlnt-l]  <  Rj) 
tap  +«■  x[anlnt]*delta_r[a 

•lse 

continue ; 


t] 


) 


for (anlnt  -  0;  anlnt  <  M;  anlnt++) 

if (vec_r [anlnt]  >■  Rj  44  Rj  <■ 
return ( log ( CKST1 ) +K* log ( x I 
else  if (vec_r [anlnt]  >-  Rj  44  Rj  >  0.35  44  Rj  <=  0.65) 
return ( log ( CKST2 ) +K*log (x ( anlnt ] ) -2 . 0*tmp) ; 

•lse  if (vec_r[ anlnt]  >■  Rj  44  Rj  >  0.65  44  Rj  <*  1.0) 


>  Rj)) 
hint] ; 


Sint] ; 

4  (Rj  <  vec_r [anlnt] ) ) 
Int]; 


0.35) 
anlnt] ) -2 . 0*tmp) ; 


) 


return ( log ( CKST3 ) +K*log (x 


/**********  ji  *******  **********  *******  ********** *************** ************/ 
/* 


anlnt])-2.0*tmp) ; 


{ 


V 

V 

V 

V 

V 

V 


/*  Name  :  euclid_norm(num_elem, array) 

/* 

/*  Purpose;  Compute  the  euclidean  nonji  of  vector  array []  of 
/*  nuaelem. 

/* 

/  ******  ******  **********  mum****  ******  **■»**  *******************************/ 

double  eucl id_norm ( num_el ea , array) 

int  num_elem; 
double  array[]; 


register  int  anlnt; 
double  norm; 


norm 


0.0; 


for (anlnt  -  0;  anlnt  <  num_elea;  anlrit++) 
norm  +«  array [anlnt] *array [anlnt] ; 


return (sqrt (norm) ) ; 


I 

j 
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